Compressibility of graphene 
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We develop a theory for the compressibility and quantum capacitance of disordered monolayer 
and bilayer graphene including the full hyperbolic band structure and band gap in the latter case. 
We include the effects of disorder in our theory, which are of particular importance at the carrier 
densities near the Dirac point. We account for this disorder statistically using two different averaging 
procedures: first via averaging over the density of carriers directly, and then via averaging in the 
density of states to produce an effective density of carriers. We also compare the results of these two 
models with experimental data, and to do this we introduce a model for inter-layer screening which 
predicts the size of the band gap between the low-energy conduction and valence bands for arbitary 
gate potentials applied to both layers of bilayer graphene. We find that both models for disorder 
give qualitatively correct results for gapless systems, but when there is a band gap in the low-energy 
band structure, the density of states averaging is incorrect and disagrees with the experimental data. 



I. INTRODUCTION 

Recently, graphene has become a highly-studied elec- 
tronic system with many potential uses in electronic de- 
vices, optics, and sensing applications^ While significant 
theoretical effort has been concentrated on the transport 
properties^ (such as the minimal conductivity, localiza- 
tion effects, and signatures of relativistic behavior of the 
charge carriers), bulk thermodynamic quantities such as 
the compressibility of the electron liquid in monolayer 
and bilayer graphene have also received substantial at- 
tention. The compressibility is an important quantity 
because it is possible to extract information about the 
electron-electron interactions and correlations directly 
from these measurements, and to gain information about 
fundamental physical quantities such as the pair corre- 
lation function. The quantum capacitance 3,4 is also an 
important quantity in the context of device design and 
fundamental understanding of the graphene system. In 
fact, the compressibility is often inferred from the mea- 
surement of quantum capacitance. 

Experiments^ have shown that the contribution to 
the compressibility of monolayer graphene from electron- 
electron interactions is small (amounting to a 10 — 15% 
renormalization of the effective Fermi velocity). This was 
explained within the Hartree-Fock approximation^ Ca- 
pacitance measurements of bilayer graphene have also 
been carried out recently*^ and they purport to show 
data which indicates that the quadratic approximation 
of the band structure of this material does not produce 
the correct predictions for the quantum capacitance of 
the graphene sheet. Other measurements which claim 
to access the quantum capacitance directly have shown 
similar behavior^ 

Theoretical work has been published which considers 
the compressibility of monolayer— and bilayer— graphene 
at the Hartree-Fock level and in the random phase 
approximation 1 ^ (RPA). They find that interactions are 
more important in the bilayer system (due to the finite 
density of states at charge- neutrality) , but that in the 



RPA the contributions from exchange and correlation 
almost cancel each other out leaving a small positive 
compressibility at all densities. A consensus has devel- 
oped that many-body corrections to the compressibility 
in graphene are at most of the order of 10% at reasonable 
carrier densities^ 

In this article we discuss theoretically the compress- 
ibility of monolayer and bilayer graphene, and specif- 
ically employ the four-band model of bilayer graphene 
(which generates the non-quadratic band structure). We 
will highlight the similarities and differences between this 
model and the linear and quadratic cases, and take into 
account the possiblility of a gap at the charge-neutrality 
point. We find that the non-quadratic band structure sig- 
nificantly affects the theoretical predictions for the com- 
pressibility and brings them into qualitative and semi- 
quantitative agreement with experimental findings. We 
also include disorder in the form of electron-hole pud- 
dles created by charged impurities in the vicinity of the 
graphene. We discuss two different procedures for includ- 
ing this effect statistically and determine the validity of 
each model by comparing with experimental data. The 
first assumes that the disorder creates a spatial variation 
in the density of carriers which can be modelled by av- 
eraging a physical quantity such as compressiblity over a 
range of densities. The second model assumes that dis- 
order affects the local density of states which can be av- 
eraged to give an effective bulk density of carriers which 
can then be used to compute the physical quantities of 
interest. We find that for gapless graphenes, inclusion of 
disorder by either method predicts qualitatively correct 
results in agreement with experiments, but when a gap is 
opened, the stage at which the averaging is done makes 
a critical difference. 

Therefore, we have dual complementary theoretical 
goals in this work. Our primary goal is to develop an 
accurate (and essentially approximation-free) theory for 
the graphene compressibility within the non-interacting 
electron model which incorporates all aspects of the re- 
alistic band structure as well as effects of the disorder 
within reasonable and physically motivated approxima- 
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tions. Our secondary goal is to compare our theory 
closely and transparently with the existing experimen- 
tal data on the compressibility of bilayer graphene so 
as to obtain some well-informed conclusions about the 
possible role of many-body exchange-correlation effects, 
which have recently been studied the the literature 1 ^— 
but which are beyond the scope of the current work. We 
believe that an accurate calculation of graphene com- 
pressibility within the free electron theory neglecting all 
interaction effects is warranted for (at least) three rea- 
sons: (?) without an accurate one-electron theory, quan- 
titative conclusions regarding many-body effects are not 
useful; (ii) many-body theory is, by definition, approx- 
imate and, as such, subject to doubts; (Hi) many-body 
corrections to graphene compressibility are claimed to be 
small with considerable cancellation between exchange 
and correlation effects^. In comparing our theoretical 
results with the experimental data, we find that the in- 
clusion of disorder effects in the theory is of qualitative 
importance, particularly at low carrier densities near the 
charge-neutral Dirac point (where the many-body effects 
are expected to be largest). The inclusion of disorder 
through two alternative physical approximations and the 
eventual validation of one of the models of disorder by 
comparing with the existing experimental data is an im- 
portant accomplishment of our theory. 

To outline the structure of this article, we begin in the 
remainder of this section by giving details of the band 
structures of the three types of graphene which we con- 
sider. In Section|TT]we introduce a model for determining 
the band gap in bilayer graphene with arbitrary gate po- 
tentials applied to both layers. Then, in Section Mil we 
describe the calculation of K = |^ in the homogeneous 
(non-disordered) situation (where /i and n are the chem- 
ical potential and carrier density respectively) for each 
of these three cases and discuss the main features of this 
quantity in relation to the quantum capacitance and com- 
pressibility. In Section IIVI we introduce the two models 
for the inclusion of disorder and present the fundamental 
theoretical results for each. Then, in Section [V] we com- 
pare the predictions of our theory to experminental data 
given by capacitance measurements of bilayer graphene 
which include the quantum capacitance of the 2D elec- 
tron gas. Finally, in Section fVTl we summarize our results 
and give some general discussion of our findings. 

In the following, we discuss three different models for 
the single-particle dispersion of graphene which we de- 
note by Ek- We shall refer to these three cases as lin- 
ear, quadratic, and hyperbolic graphene to distinguish 
them. For the monolayer, we have E k = Tivpk as is well- 
known, and we indicate all quantities which follow from 
this dispersion by attaching a superscript T as we have 
here. For bilayer graphene, there are two commonly- 
used approximations for the low-energy single-particle 
dispersion. The simplest is the quadratic approximation 
E q k = h 2 Vpk 2 /^/i = h 2 k 2 /(2m) where 71 is the inter- 
layer hopping parameter from the tight-binding theory 
and m = 7i/(2vp) is the effective mass of the electrons. 
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FIG. 1. Sketch of the 'sombrero' dispersion in hyperbolic 
graphene with u > 0. The depth of the minimum is ex- 
agerrated for the purposes of illustration. The various wave 
vectors and energies that are labeled are defined in the text. 



This dispersion comes from a low-energy effective theory 
of bilayer graphene 1 -^ which essentially discards the two 
split bands and confines electrons to those lattice sites 
not involved in the inter-layer coupling. This quadratic 
bilayer dispersion is assumed to hold for low carrier den- 
sities, i.e. at low Fermi energy 3 Alternatively, the full 
tight-binding theory 14 yields the single-particle disper- 
sion 




h 2 v 2 F k 2 



- + h 2 v 2 F k 2 ( 1 2 + u 2 ) 

(1) 

for the low-energy bands, which is illustrated for the con- 
duction band in Fig. [TJ This dispersion allows for the 
inclusion of a band gap parameterized by the energy u 
which corresponds to the potential energy difference be- 
tween the two layers. The 'sombrero' shape of the dis- 
persion implies that the minimum gap is found at wave 

vector k* — 2 ™ vp an d the energy of the mini- 

mum is therefore E* = — 7 1 " This non-monotonic 

dispersion relation leads to some important features in 
the Fermi surface at low density. For medium and high 
density regimes (when kp > ) , the Fermi wave vector 

is given by &f = V 7m an d the Fermi surface is disc- 
shaped. However, if y/irn < (as sketched in Fig. [T]) 
then the Fermi energy is in the 'sombrero' region and the 
Fermi surface is ring-shaped 15 with kp- < k < kp + and 
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7x + u 2 



u 2 ±2TTh 2 v 2 F n. (2) 



Clearly these characteristic features of the full hyper- 
bolic bilayer band dispersion are lost in the linear ( "high- 
energy") and quadratic ("low-energy") approximations 
often adopted in the literature for simplicity. 

The parameter u can be treated in two ways. Initially 
(and in Section |IV[) . we consider it to be a phenomeno- 
logical parameter which can be chosen arbitrarily to rep- 
resent a finite gap at charge-neutrality. However, in or- 
der to accurately compare our theory with experiment 
we must determine a method of relating u to the exter- 
nal gate potentials. To do this accurately, the screening 
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of the external field by the charge on the two layers of 
bilayer graphene must be taken into account, and we de- 
scribe this process in Section HT1 



II. INTER-LAYER SCREENING 

In order to make an accurate comparison with exper- 
imental data, we must properly map the gate voltages 
applied in experiment to the parameters n and u which 
appear in our theory. We assume the situation sketched 
in Fig. [2ja). The density is easily computed via elemen- 
tary considerations to be n — eo (VbCb / db + VtCt / dt) / e 
where V is the voltage applied to a gate, d is the width 
of the gate layer, e is the dielectric constant of the gate, 
and the subscript b and t respectively denote the bottom 
and top gates. The zeroeth-order approximation for the 
on-site potential u is to say that the electric field cre- 
ated by the gates directly gives the electric field E in the 
graphene. The potential is then Utj> — ieEd/2. How- 
ever, this approximation does not take into account the 
screening effect of the charge on the two graphene lay- 
ers. Previously, various authors have investigated this 
problem and shown that screening reduces the size of the 
gap 14 i 16 ' 17 . We extend this work to the case where the 
total density is finite and gates are present both above 
and below the graphene layer. 

Since the screening depends only on the electric field 
created by the gates and the charge densities on the two 
layers, it is independent of the disorder caused by charged 
impurities. Both the screening and the disorder effects 
are determined by the net carrier density, and both leave 
this quantity unchanged. Therefore, we can use the fol- 
lowing analysis to set the size of the gap for any given set 
of external conditions, and use the resulting band struc- 
ture to calculate the effect of disorder for that situation. 
This is allowed because neither disorder nor the screening 
induces any extra charge in the graphene. Additionally, 
we will assume in Section |IV] that the disorder potential 
is the same in both layers meaning that it does not in- 
duce any inter-layer scattering and therefore leaves the 
net carrier density in each layer unchanged. 

According to Gauss' law, the electric field directed 
away from an infinite-sized, positively charged sheet is 
E = qn/2eQ, where q — — |e| is the charge of the carriers 
on the sheet. Hence, we can write the internal field (i.e. 
the field between the graphene layers) as 



e 

E = E cxt + —(ra t - n b ) 
Zeo 



(3) 



That this internal field is different from the external one 
causes the free charge on the two layers to rearrange, 
which in turn further modifies the internal field. A self- 
consistent procedure can be constructed to find the final 
solution for the layer charge densities and the internal 
electric field. This additional field causes an additional 
Hartree energy of ut — Ub — edE to be added to the on- 
site potential in the original Hamiltonian. This modified 
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FIG. 2. (a) The geometry of the screening problem. The 
external field i5 cxt = (Eb + Et)/2 is screened by redistribution 
of the charge density n = nb+nt between the two layers of the 
graphene. (b) The screened gap u as a function of the bare gap 
Wext for various densities, (c) The screened gap as a function 
of the top gate voltage for various fixed back gate voltages 
and the experimental parameters used in the experiments of 
Young et al. 7 '^ s Back gate voltages are quoted in Volts. 



on-site potential will be smaller than in the unscreened 
case, and therefore the gap will also be smaller. 

We implement this procedure for a given pair of gate 
voltages Vb and V t as follows. The first step is to compute 
the unscreened on-site potential in the two layers, which 
is given by the expression above with n t = n b so that E = 
E cxt . This gives u eKt = edE cxt = ed(V b /d b - V t /d t )/2. 
We assume that the potential is symmetrical so that 
Ut = u/2 and Ub = —u/2. This determines the single- 
particle Hamiltonian, from which we can compute the 
wave functions. The layer density is then found by sum- 
ming the wave function weight in each layer: 
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states 



(10 



AL 



M 2 ) 



where A and B label the two inequivalent lattice sites in 
each layer. The sum runs over all valley, spin and wave 
vector states. A similar expression is given for the lower 
layer by substituting t — > b throughout. Computing these 
densities gives all the information required to evaluate 
the internal electric field, from which we can find the 
modified on-site potentials via Eq. ([3|). These potentials 
are then used to calculate the modified layer densities, 
and this cycle is repeated until convergence is reached. 
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Figure (2jb) shows the screened gap as a function of the 
external gap for various density of carriers. The density 
plays a role because the more free carriers there are, the 
more effectively the external field can be screened and the 
smaller the screened gap is. This is shown in the figure 
because the magnitude of the gap for \n\ = 5 X 10 12 cm -2 
is always smaller than that for \n\ = lx 10 12 cm~ 2 . Figure 
[He) shows the screened gap as a function of top gate volt- 
age assuming that the back gate voltage is fixed. We have 
used the parameters from the experiments by Young et 
alZ*^. The peak in each trace is caused by the approach 
of the system to charge neutrality where there are no 
excess carriers and hence the external electric field is un- 
screened. At this point, the screened gap is equal to the 
bare gap. When there is a finite density of excess car- 
riers, the external field is screened and the gap reduced 
accordingly. As the back gate voltage is changed, the 
position of the charge neutrality point shifts and the size 
of the gap at charge neutrality increases. Inclusion of 
inter-layer screening turns out to be important for un- 
derstanding the experimental data. 



III. 



CALCULATION OF ^ 



an integral over k. The indefinite integral in question is 
given by 



2A 



kEl dk. 



(5) 



To find the total energy, we write the sum over states dis- 
cussed above as an integral over the wave vector. Then, 
evaluation of the Fermi function naturally leads to an 
expression containing I(k) at the limiting values of the 
wave vector for the specific value of the density: 

2A f°° 

£ = — kE£n F (k)dk = I(k = k F+ )-I(k = k F _). 
n Jo 

The chemical potential can be calculated by transforming 
the derivative into one for the limiting fc: 



^ A dn 



1 d£ 1 dk F+ dI(k F+ ) 



A dn dk. 



1 dk F _ dI(k F _ ) 
A dn dkp 



However, it is clear from the fundamental theorem of 
calculus and the definition of / that dl/dk = 2AkE k /n so 
that the second derivative is straightforward to calculate 
and we have 



In this section, we describe the calculation of & for 
the homogeneous (non-disordered) system. The results 
derived here will be used later to compute the same 
quantity for a disordered system. The single-particle dis- 
persions Ek of monolayer and bilayer graphene are well 
known (as discussed in Section UJ. From these relations, 
we can compute the total kinetic energy for excess charge 
carriers by summing the energy of carriers in filled states: 
E = Estates -^fe n F(^) where n F {k) is the Fermi distribu- 
tion function with wave vector k — |k|. Then, the chem- 
ical potential is defined to be the change in energy with 
the addition of a particle, which is expressed as a = \^ 
where n is the carrier density and A is the area of the 
graphene flake. From this we can compute K = ^p, and 
the quantum capacitance Cq and compressibility k are 
then linked to this quantity via 
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, dn 
da ' 



dn 



(4) 



The linear and quadratic band structures lead straight- 
forwardly to the following expressions 



K> 



and K q 



h 2 V F TT 
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The linear dispersion leads to a square-root divergence at 
zero density while the quadratic version gives a constant 
for all densities. It is possible to obtain an analytical 
expression for K h in various limits and we present them 
in the Appendix. Here, we describe a process for writing 
down the answer. Starting from the expression for the 
total energy, the sum is evaluated by transforming to 




(6) 



The remaining derivatives are simple to compute from 
Eqs. (Q} and ([2]). As described previously, when \E F \ < 
u/2 the values of k± are given by Eq. $Z§. When \E F \ > 
u/2 (which is trivially the case for u = 0) then k- = 
and k+ = \/im and only the first two terms contribute to 
K h . This is the central result of this section, but in the 
Appendix we have evaluated this expression for various 
limiting cases. 

The three versions of K are plotted in Figure [3ta) 
where the density is assumed to be a tunable external pa- 
rameter and the gap is fixed. In the limit of small n, the 
gapless hyperbolic band structure is roughly quadratic 
and so K h is very close to K q . However, the linearity 
of the hyperbolic dispersion quickly asserts itself to re- 
duce the size of K h until, at high density, it approaches 
the value of linear graphene. Thus K h in the ungapped 
regime interpolates between these two limiting cases as 
one expects. We also show K h for the gapped case (red 
or light gray line) where we take u = 50meV. When the 
density is high and the Fermi energy is well above the gap 
region, K h tends to the same functional form as for the 
ungapped scenario. However, when the density is low, 
K h goes to zero for any finite u. The large reduction in 
K h by the gap will be an important feature for the com- 
parison with experiment. In Figure [3jb) we show the 
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FIG. 3. (a) The inverse density of states, K = for 
linear, quadratic, gapless hyperbolic and gapped hyperbolic 
(it = 50meV) dispersions, (b) K for different gap sizes at low 
carrier density. 



K h at low carrier density for several different values of 
the gap. One immediately notices a discontinuity when 
n = u/(hvp) which is caused by the change in topology 
of the Fermi surface from a ring to a disc. When a gap 
is present, there is also a (5-function divergence at n = 
(not shown in the plots) caused by the step in the chem- 
ical potential as the Fermi energy goes from the valence 
band to the conduction band. The prefactor of this di- 
vergence is the size of the gap so that it is not present 
when u = 0. 



IV. MODELS FOR DISORDER 

We now discuss the effects of disorder caused by 
electron-hole puddles which are formed by charged im- 
purities in the environment We will describe two in- 
equivalent methods of averaging over disorder and high- 
light the differences between them. In what follows, we 
must be careful to distinguish global (averaged) quanti- 
ties from their local counterparts. Therefore, we use an 
overbar to denote a quantity which has been averaged 
over disorder and which is therefore global. We start by 
assuming that charged impurities (which may be located 
in the substrate or near the graphene) create a local elec- 
trostatic potential which fluctuates randomly across the 
surface of the graphene sheet. In the case of the bilayer, 
we assume that this potential is felt equally by both lay- 
ers. We then make the transition from the spatial varia- 
tion of the potential to an average by assuming that the 
value of this potential at any given point can be described 
by a statistical distribution. 20 This distribution can then 
be used to average some pertinent quantity to obtain the 
bulk result for the disordered system. But it is not clear 
a priori what the optimum method of performing this 
average is. In order to discuss this question, we use two 
different methods and (in the next section) compare the 
results of each to experimental data. 

The first model assumes that the disorder potential 
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FIG. 4. Electron, hole and total carrier densities for linear, 
quadratic and hyperbolic graphene for s — 20meV. 



directly affects the density of carriers in the graphene so 
that the density is distributed according to a Gaussian 
distribution P parameterized by width v 



P(n) 



1 



exp 



n 

2^2 



Then, the overall value of K = ¥ can be found by 
evaluating the convolution of the homogeneous (non- 
disordered) K with this distribution: 



dfx 
dn 



dn(n — n') „, ,. , , 

-JT K-P(n')dn'. 

a(n — n ) 



(7) 



This has the effect of broadening any sharp features in 
the homogenous K, such as the 5-function at n = in 
a gapped system or the step associated with the change 
in topology of the Fermi surface in the hyperbolic band. 
Specifically, evaluating this convolution gives the follow- 
ing contribution to K for the hyperbolic band from the 
(5-function: 

U71 exp (-n 2 /2(Sn) 2 ) 



sju 2 + a 



2tt 



so that this spike is broadened into a peak where the 
height is determined by the size of the effective gap 2E* . 
This model corresponds to that taken in Ref. (a, except 
that our theory includes the inter-layer screening which 
reduces the gap size. 

The second model assumes that the disorder potential 
introduces variation in the density of states (DoS). The 
potential itself is assumed to follow a Gaussian distribu- 
tion P(V) parameterized by width s, and therefore the 
disordered DoS is given by^ 



D(E) 



D (E -V)P(V)dV 



(8) 



where Dq(E) is the homogeneous (non-disordered) DoS. 
This modified DoS then gives the zero-temperature elec- 
tron density by integrating over energy so that 



/oo 
D(E)&(E F - E)dE 
-OO 



(9) 
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where Ep is the global Fermi energy and Q(x) is the 
Heaviside step function. A similar expression exists for 
the density of holes. This density is assumed to hold 
throughout the graphcne so that the averaged density can 



now be used to evaluate K 



~T7 _ dp{n) 



dn 



Since we assume that 



the disorder does not induce any additional charge in the 



2D system, the net carrier density n e 



is conserved. 



Therefore, to compute the effective Fermi energy Ep in 
the presence of disorder we solve the equation rf(Ep) — 
n h (E^) = nl(E F ) - n%(E F ) for ~Ep~. However, the total 
carrier density n e + n h is not conserved, and the primary 
effect of disorder in this model is to introduce a finite 
total density for any finite amount of charged impurities. 

The integrals in Eqs. (|8]) and ([9]) can be computed an- 
alytically for the linear and quadratic bands if we assume 
that P(V) is a Gaussian distribution. They are 



and 
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2nh 2 vl 



± E F erfc(Ty) 



where the upper sign corresponds to the electron density 
and the lower sign to the hole density, erfc is the com- 
plementary error function, and y = Ep j \\[2s) with s the 
standard deviation of the distribution of the disorder po- 
tential. We see that the effect of disorder is to introduce 
a tail in the electron density which extends into the con- 
duction band, the size of which increases with increasing 
disorder strength. In the hyperbolic case, the integrals 
in Eqs. © and © cannot be evaluated analytically, so a 
numerical comparison must be done with the other two 
cases. 
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erfc(=Fy) ± -^=e y 
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The expressions which must be computed for the hyperbolic band are 
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s\[2~k irfi 2 v 2 7 



dE 



E-E* 



dV(E - V) [2 + T(E - V)] 



,-V 2 /(2s 2 ) 



E-E* 



dV (E -V)[2- T{E - V)] e 



and 



n hh 



1 I 

s\[2~k Trfh 2 v F Je f 



poo r />oo 

/ dE - dV{E -V)[2 + T{E - V)} e 

J Ep L " E-\-E* 



E+i 



(E - V) [2 - T(E - V)] 



,-V 2 /(2s 2 ) 



E+E' 



where A = I + u 2 /7i\ ^( x ) = Ay 'y ; Ax 2 ' — 8 2 and in 
both cases the second term comes from the inner Fermi 
surface for E* < \E\ < u/2. 

Figure [5] shows the density of carriers in these three 
cases for a modest value of disorder (s = 20meV) and 
also for u = 50mcV in hyperbolic graphene. The small 
density of states in linear graphene gives a small carrier 
density at low Fermi energy. Quadratic and hyperbolic 
graphene have finite density of states at charge-neutrality 
so the density is higher in this case. The density in hyper- 
bolic graphene increases slightly quicker than quadratic 
graphene because of the linear increase in the density 
of states compared to the constant density of states in 
quadratic graphene. Also, the finite density of states at 
the charge-neutrality point in quadratic and hyperbolic 
graphene means that the effect of disorder is significantly 
stronger in these systems than in linear graphene. In fact, 
substituting Ep = into the total density gives the de- 
pendence of the density at the charge-neutrality point as 



a function of the disorder strength s as 

so that the linear increase in n q (Ep = 0) is faster than 
the quadratic one in n l (Ep = 0) for small s. The density 
at charge-neutrality in hyperbolic graphene is identical to 
that in quadratic graphene because the density of states 
is the same in both cases. We emphasize that both mod- 
els of disorder considered here are physically motivated 
and a priori there does not appear to be any practical 
reason to choose one over the other. In the next section, 
we compare theory with recent graphene compressibility 
experimental data to establish the comparative validity 
of our disorder models. 
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FIG. 5. Classical circuit diagrams for the experiments by (a) 
Henriksen et aZ.— , and (b) Young et all 



EXPERIMENTAL COMPARISON: 
EFFECT OF DISORDER 



THE 



We now discuss our results in the context of recently 
published experiments. We restrict ourselves to compari- 
son with the bilayer because it is known that the effect of 
disorder in monolayer graphene is small, and this has al- 
ready been discussed in the literature. 6 In particular, two 
measurements have been reported where the capacitance 
of bilayer graphene has been measured and the compress- 
ibility extracted from those measurements The quan- 
tum capacitance is linked to ^ since Cq — Ae 2 where 
A is the surface area of the sample. The two experiments 
actually measure slightly different capactitances which 
are represented by the two circuit diagrams in Fig. [5] and 
correspond to the following expressions for the measured 
capacitances Ch and Cy. 



Cy= m^ +Cs 



e t d„ + dt 



C 



H 



e e b e t d g A 



dg(e b d t + e t db) + d b d t 



Young Ref . 7 
Cs Henriksen Ref. 



The notation follows that defined in Fig. HJa), except 
that d g is an effective length associated with the quantum 
capacitance as 



C Q = 



Ae 



dg = 



fo 

e 2 dn 



The quantity Cs is the background (or stray) capaci- 
tance associated with, for example, the substrate and 
experimental equipment. For relevant experimental pa- 
rameters, d g <C db, dt so that Ch oc d g cx We find the 
relevant experimental parameters from the papers^ and 
use them to produce the plots shown in Figs. [5] and [7J 
Note that if the quadratic band structure is used to com- 
pute Ch and CY, so that d g = h 2 Vp€ n/(e 2 ji) then there 
is no density dependence in either Ch or Cy, which is 
qualitatively different to the experimental data in Figs. [6] 
and [7J Therefore the full band structure is an essential 
part of our theory. 

Figure [S] shows the comparison of data from both dis- 
order models to the experiments^ of Young et al. In this 
case, the back gate is kept at a fixed potential while the 
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FIG. 6. Comparison of theoretical results including disorder 
for the measured capacitance Cy with experimental data by 
Young et a/.— We have taken 71 = 0.35eV, and vp = 1.05 x 
10 6 ms -1 . 



top gate is varied. This has the effect of simultaneously 
changing the size of the gap and the carrier density so in 
this case we treat these two quantities seperately and use 
the full theory of the inter-layer screening described in 
Section [XT] to compute the size of the effective gap. When 
the gap is small at low density {i.e. for V b = 10V, plot 
(a)) the both theories of disorder predict qualitatively 
correct answers. The difference between the experimen- 
tal and theoretical curves can be reduced by changing 
the value of the stray capacitance as a fitting parameter. 
For very low values of disorder, a small spike is present 
at zero density in the density averaged data. This is due 
to the very small band gap that is present and the asso- 
ciated (5-function in . Larger values of disorder smear 
this divergence out into a dip, the shape of which matches 
the experimental data well. The DoS averaged plots do 
not show this dip because carrier density is always fi- 
nite in this model, and the effect of increasing disorder is 
to increase the minimum value of the carrier density thus 
decreasing the depth of the dip. When the gap is large at 
zero density (i.e. for large back gate voltage, Vb = 70V, 
plot (b) and V b = -80V, plot(c)) then the DoS aver- 
aging procedure predicts qualitatively wrong results for 
low values of disorder. There are two main features in 
the low density part of these plots. First, for finite gap, 
Fig. [3] shows that ^ tends to zero as density decreases. 
This gives a peak in the measured capacitance (because 
Cy cx 1/(|£ + f^-)) which is blurred by finite disor- 
der. The second feature is the (5-function at the band 
edge (or, equivalently, at zero density) which is broad- 
ened by finite disorder in the density averaging model, 
but which plays no role in the DoS averaging model be- 
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FIG. 7. Comparison of theoretical results including disorder 
with experimental data by Henriksen et oZ.— (a) u cxt = 
compared with the density averaged model for zero gap and 
a small gap; and u — 25meV for v = 3x 10 11 cm -2 . (b) it ex t 
finite, compared with the density averaged theory for several 
values of the gap. We have taken v — 3 x 10 11 cm~ 2 and the 
legend is as shown in (d). (c) DoS averaged theory compared 
with it cxt = expermimental data for s — and s — 30meV. 
(d) Mext finite compared with DoS averaged theory for two 
values of the gap. s — 30meV in each case. We have taken 
71 = 0.35eV, and vf = 1.05 x 10 6 ms _1 throughout. 



cause this model always gives rise to a finite density of 
carriers. This means that the DoS average cannot pre- 
dict correct results in the low-density, finite gap regime. 
However, the density averaging model does include the 
effects of the ^-function spike by sampling many values 
of density, and crucially including n = 0. This spike in 
^ manifests as the dip in the measured capacitance and 
parameters can be chosen so that the size and shape of 
the dip match the experimental data quite closely. We 
therefore believe the density-averaging procedure of Eq. 
([7|) to be the preferable model for including disorder over 
the the DoS-averaging model. 

This general picture is also seen in the comparison of 
the two disorder theories with data by Henriksen et al. 
as shown in Fig. [7] Panels (c) and (d) show the compari- 
son of the experimental data to the DoS averaged model. 
In this experiment, the measured capacitance is propor- 
tional to , so in the gapped case in panel (d) , the step 

in due to the change in topology of the Fermi surface 
is clearly seen for large gap. The zero gap data in panel 
(c) shows qualitatively the correct features, although the 
stray capacitance must be altered to get a better fit and 
the peak is taller in the experimental data than in the the- 
ory. But when the gap is opened as shown in part (d), the 
DoS averaged theory predicts qualitatively wrong results 



because it does not account for the broadened (S-function 
associated with the band edges. The experimental data 
also shows a distinct asymmetry between the electron 
and hole density regimes. This is not explained in our 
model, but can be added by allowing for an asymmetry in 
the underlying band structure, for example by including 
specific next-nearest neighbor hops in the tight binding 
model. 21 On the other hand, the density-averaging the- 
ory (shown in panels (a) and (b)) do predict the correct 
behavior. For the zero-gap experiment, the variation of 
the predicted Cr is too small, and better agreement can 
be found by allowing for a small gap to open as shown by 
the dashed line. In panel (b), it is clear that a gap size of 
sa 40meV gives a semi-quantitative fit to the experiment. 



VI. CONCLUSION AND DISCUSSION 

In conclusion, we have presented a full calculation of ^ 
(and hence the quantum capacitance and compressibil- 
ity) for graphene and have discussed two different meth- 
ods for including disorder effects. These models require 
averaging the effect of the disorder potential in either 
the density of states or directly as a variation in the local 
electron density. We find that when there is no gap in 
the system, these two models predict qualitatively similar 
results. But when there is a gap present, the DoS averag- 
ing method fails to account for the step in the chemical 
potential due to the band edges and therefore predicts 
qualitatively incorrect results. So, we present two conclu- 
sions: First, that the single- particle theory for the com- 
pressibility of graphenes appears to allow qualitatively 
good predictions for thermodynamic quantities such as 
the compressibility. Second, the inclusion of disorder via 
statistical averaging must be done carefully because the 
choice of where the averaging is done critically affects the 
outcome. We find that the averaging should be carried 
out at the level of the experimental quantity rather than 
the density of states. 

We believe that the effects left out of the theory, par- 
ticularly the many-body exchange-correlation corrections 
are likely to be qualitatively unimportant for under- 
standing the experimental data for currently available 
graphene samples which are dominated by disorder. In 
particular, disorder effects are strongest at the lowest car- 
rier densities where exchange-correlation effects become 
more important quantitatively. This makes an observa- 
tion of many-body effects through compressibility mea- 
surements quite challenging. 

We thank J. P. Eisenstein, E. A. Henriksen, P. Kim, 
and A. Young for discussions of their experimental re- 
sults. This work was supported by US-ONR and NRI- 
SWAN. 



9 



Appendix: Analytical results for K h 



When u > 0, \E F \> u/2 



In this Appendix, we give the derivation of K h , and 
write down the full expression for it in various cases. The 
full form of the indefinite integral in Eq. is 



m = 2 AJ kE H 
All 



dk 



1 



irh 2 v 2 2A^A 



7 2A (8x-<Z>- 



+ 12S 2 log 2A (-A 

where x = (hv Fk j '71) 2 , 6 = u/(2^i), 

and A = 1 - 



1 - AS 2 ) 3 
+ V2AE 



$ = 
-AS 2 . 



Explicitly 



S = V2a; - $ + 1 + 26 
differentiating with respect to k verifies the Fundamental 
Theorem of Calculus which we utilized in the description 
of Eq. ([HI). In the general case, the derivative of the 
energy with respect to the wave vector is 



When \Ep\ > u/2 then k F — \prni and k F _ = 
so that the derivatives with respect to n are elementary 
and only the first two terms of Eq. ([6]) contribute to K h . 
Therefore 



K h {\E F \ > 



7i 



2n Q 



471 



1 j n_ tiGu 

4 ' n n 



1 - 



n 0u /n a 



X 1 _n_ np u 

4 ' n n 



where n 0u = (7^ + u 2 ) / (irfi 2 v F ) . 

When u > 0, \E F \ < u/2 

However, when \E F \ < u/2 we need to use the expres- 
sions for the two Fermi surfaces in Eq. © so that 



dE h kF 
dk F 



h 1 v 2 F k F 
E L 



7? + u 2 



n 2 v 2 F k 2 F ( 1 2 + u ^) i 



(A.l) 



The u = limit 



and hence 
d 2 k 



dn 2 



dk 



F± 



dn 



ik 



Ak 



1 

n 0u 



n 
n 0u 



± 1 



Ak 2 
^F± 



± 1 



n u 



im for all densities. Therefore the 



In this case, k F 
derivatives of k F with respect to n are elementary, the 
k F terms in Eq. (jB]) trivially do not contribute and we 
have 



K h — 



7i 



"0 



Substituting these expressions for the derivatives into Eq. 
© gives the simplified expression 



K h (\E F \ < u/2) 



1 / n 








_k F \no u 





n 

n 0u 



2 dE u 



dkr 



2 dE u 



dkr 



where n = "/ 2 /(nh 2 vj 



where the derivatives of the Fermi energy with respect to 
the Fermi wave vector are given by making the appropri- 
ate substitutions in Eq. (|A.1|) . 
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